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ABSTRACT 

A  bi-quadratic  isoparametric  plate/shell  bending  finite  element  is  developed 
to  study  the  behavior  of  isotropic  and  laminated  composite  plates.  The  element 
is  based  on  Mindlin-Reissner's  theory  and  the  principle  of  virtual  displacements. 
The  element  is  implemented  in  a  computer  program.  Results  are  presented  and 
compared  with  analytical  solutions  to  validate  this  element.  Good  agreement  is 
observed  for  thin  plates,  while  discrepancies  are  noted  for  thick  plates.  Effects  of 
various  integration  schemes  on  the  element  performance  are  presented.  Convergence 
studies  for  laminated  composites  for  different  fiber  orientations  are  also  discussed. 
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I.  INTRODUCTION 

A.     THE  SCOPE  OF  THE  THESIS 

Finite  element  analysis  provides  a  general  tool  to  solve  problems  in  struc- 
tural mechanics.  The  methodology  is  applicable  for  static  and  dynamic  response  of 
structures  and  in  predicting  the  elastic  stability  limits. 

The  focus  of  the  present  study  is  to  develop  tools  to  analyze  laminated  com- 
posite plates  and  validate  the  model  by  comparing  with  known  solutions. 

More  specifically,  the  objectives  of  the  present  study  are: 

a)  to  review  some  of  the  pertinent  literature  in  the  area  of  laminated 
composite  plates. 

b)  develop  a  finite  element  for  the  analysis  of  composite  plates. 

c)  develop  consistent  mass  and  load  matrices. 

d)  study  the  effect  of  thickness  to  characteristic  length  of  the  plate. 

e)  study  the  effects  of  integration  schemes. 

The  outline  of  the  remainder  of  this  thesis  is  as  follows: 

The  basic  formulation  of  the  stiffness,  mass  and  load  matrix  for  the  bending 
of  flat  plates  using  laminated  composite  materials  are  described  in  Chapter  II. 

Chapter  III  addresses  certain  aspects  of  computational  implementation. 

Chapter  IV  describes  some  test  cases,  example  calculations  and  comparison 
with  classical  plate  theory. 

Finally,  Chapter  V  reflects  experience  gained  and  some  suggestions  for  future 
research. 


B.     LITERATURE  SURVEY 

The  finite  element  method,  [Ref.  1],  may  be  described  as  a  general  discretiza- 
tion procedure  of  continuum  problems,  posed  by  mathematically  defined  statements 
with  applications  to  several  engineering  analysis  problems. 

A  brief  literature  review  pertaining  to  the  analysis  of  plates/shells  using  finite 
element  approximation  is  presented  in  the  following  paragraphs.  A  significant  con- 
tribution to  include  shear  is  given  by  Mindlin  [Ref.  2],  while  Hughes  etal,  [Ref.  3] 
adapt  this  theory  to  develop  finite  elements  for  the  analysis  of  isotropic  plates  [Refs. 
3,  4].  Laminated  plate  theory  based  on  the  classical  KirchofF  hypothesis  has  been 
established  by  Reissner  and  Starsky  [Ref.  5]  and  Whitney  and  Leissa  [Ref.  6].  The 
effect  of  reduced  integration  in  isoparametric  elements  was  presented  by  Zienkiewicz 
etal  [Ref.  7]  and  Hughes  etal  [Ref.  3]. 

The  finite  element  method  of  analysis  for  the  plate  bending  problem  including 
shear  deformation  has  been  presented  by  Pryor  and  Barker  [Ref.  8].  Mawenya 
describes  formulations  for  multi-layer  plates  [Refs.  9,  10].  The  higher  order  shear 
deformation  theory  of  laminated  composite  plates  was  developed  by  Krishna  Murty 
[Ref.  11],  and  Lo  et.al.  [Refs.  12,  13]  present  a  higher  order,  three-dimensional 
theory. 

Burt  [Ref.  14]  presented  a  higher  order  theory  and  compared  with  Pagano's 
elasticity-theory  solution  for  the  case  of  cylindrical  bending  and  a  symmetric  cross- 
ply  laminate  consisting  of  three  equal-thickness  layers.  Bending  of  simply  supported 
thick  rectangular  plates  was  presented  by  Srinivas  and  Rao  [Ref.  15].  Exact  elas- 
ticity solutions  for  some  particular  plate  bending  problems  have  been  obtained  by 
Pagano  [Refs.  16,  17,  18,  19]  and  Srinivas  and  Rao  [Ref.  20]. 


Application  of  classical  shell  theory,  including  transverse  shear  deformation  is 
presented  by  Vinson  and  Chou  [Ref.  21].  Naschie  [Ref.  22]  studied  large  deflec- 
tion behavior  of  orthotropic  composite  materials.  Srinivas  [Ref.  23]  developed  a 
refined  approximate  theory  for  the  static  and  dynamic  analysis  of  finite,  laminated, 
composite,  circular  cylindrical  shells  with  general  boundary  conditions. 

Plate  theories,  which  include  shear  deformation  has  been  developed  by  Whit- 
ney [Ref.  24]  and  Mau  [Ref.  25].  The  first  such  theory  for  laminated  isotropic 
plates  is  due  to  Yang,  Norris  and  Stravinsky  [Ref.  26].  Reddy  [Ref.  27]  developed 
a  higher  order  shear  deformation  theory  of  laminated  composite  plates.  Reddy  and 
Sandidge  [Ref.  28]  presented  mixed  finite  element  models  of  the  classical  and  shear 
deformation  theory.  The  effect  of  transverse  shear  deformation  on  bending  of  elastic 
symmetric  laminated  composite  plate  undergoing  large  deformation  is  presented  by 
Gorji  [Ref.  29]. 

Based  on  anisotropic  elasticity,  Hearmon  [Ref.  30]  and  Lekhnitskii  [Ref.  31] 
present  general  theories  for  laminates.  The  covariant  form  for  the  transformed  lam- 
ina stiffnesses  has  been  given  by  Tsai  and  Pagano  [Ref.  33].  Gibert  and  Schneider 
[Ref.  34]  directed  their  study  in  this  direction.  Noor  and  Mathers  [Refs.  35,  36]  pre- 
sented the  effects  of  shear  deformation  and  anisotropy  on  the  response  of  laminated 
anisotropic  plates. 

Nelson  and  Lorch  [Ref.  40]  compare  the  accuracy  of  various  plate  models  to 
predict  the  behavior  of  laminated  orthotropic  plates. 


II.  THEORETICAL  FORMULATION 

A.  INTRODUCTION 

In  this  chapter,  the  derivation  of  plate  finite  elements  based  on  Mindlin's 
theory  is  described.  The  principle  of  virtual  displacements  is  invoked  to  obtain 
equilibrium  relations. 

B.  THE  PRINCIPLE  OF  VIRTUAL  WORK 

In  this  section,  we  prove  that  total  internal  virtual  work  is  equal  to  total  ex- 
ternal virtual  work  and  equivalence  of  this  principle  to  the  minimum  total  potential 
energy  principle. 

In  general,  the  total  potential  energy  of  a  structural  system  is  equal  to  the 
sum  of  strain  energy  and  potential  energy. 

UP  =  U  +  V  (2.1) 

where, 


n„     :     Total  Potential  Energy  of  Structural  System 
Total  Strain  Energy 
Total  Potential  of  External  Loads 


U 
V 


The  total  minimum  potential  energy  requires  that  first  variation  of  total  po- 
tential energy  be  zero,  or 

6UP  =  0  (2.2) 

or, 

6U  +  SV  =  0  (2.3) 


in  other  words, 

SU  =  -6V  (2.4) 

This  may  also  be  written  in  a  different  form,  recognizing  U  as  being  the  work 
done  by  internal  forces  and  that  work  done  by  the  external  forces  being  equal  to  the 
negative  of  the  total  potential  energy  of  the  external  loads.  That  is, 

6Wtnt  =  6Wext  (2.5) 

which  is  a  statement  of  the  principle  of  virtual  work,  stating  that,  if  the  body  is  in 
equilibrium,  the  total  virtual  work  done  by  the  internal  forces  is  equal  to  the  total 
virtual  work  done  by  external  forces  for  arbitrary,  kinematically  admissible  virtual 
displacements. 

It  may  be  noted  that  the  form  in  Equation  (2.4)  is  restricted  to  conservative 
loadings  while  the  form  in  Equation  (2.5)  is  applicable  for  any  loading  form. 

The  total  internal  virtual  work  may  be  written  as 

SWint  =  I   {*}T{6e}d(A)  (2.6) 

J  vol 

where, 

{a}      :     Vector  of  Stresses 
{Se}     :     Vector  of  Virtual  Strains 
By  using  generalized  Hooke's  law  for  material  constitutive  relations,  stresses 

may  be  expressed  as 

M  =  lQ]{e}  (2.7) 

where  the  matrix  [Q]  contains  the  material  stiffness  coefficients. 

If  the  thickness  t  is  constant,  then  total  internal  virtual  work  takes  the  form, 


SWtnt  =   (  t{e}T  [Q]  {Se}d(vol)  (2.8) 

J  A 

In  order  to  derive  the  element  matrices,  the  principle  of  virtual  work,  which  is 
equally  applicable  to  the  element  as  well  as  to  the  total  structure,  is  applied  to  the 
element.  The  virtual  work  is  additive  and  results  in  the  virtual  work  of  the  entire 
structure  under  consideration. 

The  linear  strain-displacement  relations  are  given  by 

{e}  =  [B]{u}  (2.9) 

while  the  virtual  strains  are  given  by 

{6e}  =  [B]  {Su}  (2.10) 

The  operator  matrix,  [B],  is  dependent  on  the  shape  functions  and  their  deriva- 
tives, and  {u}  and  {Su}  are  vectors  of  displacements  and  virtual  displacements, 
respectively  of  the  element. 

On  substituting  Equation  (2.9),  (2.10)  in  Equation  (2.8),  we  obtain, 

SWint  =  f  t({Su}T[B}T)[Q]([B){Su})dA  (2.11) 

J  A 

or,  rearranging 

6Wtnt  =  {6uf(t  f  [B]T[Q)[B]dA)  (2.12) 

J  A 

6Wint  =  {6u}T[K}  (2.13) 

where  [K]  =  t  fA[B  ][Q][B]dA  is  the  element  stiffness  matrix. 


In  order  to  derive  mass  and  load  matrices,  we  start  from  the  expression  for 
external  virtual  work, 

SWext  =  I  T{6u}ds  -  J  x{Su}d(vol)  +  {Su}T{F}  (2.14) 

where,  T  is  the  external  force  per  unit  length  along  the  boundary  of  the  element, 
x  is  the  body  force  per  unit  volume  (inertia  force) 

{F}     is  the  point  loads  applied  at  nodal  points. 
On  substituting  for  displacements  in  terms  of  nodal  displacements  using  shape 

functions,  we  obtain, 


6Wext  =  {6u}T(  f[N]Tds)  -  {Su}T(  I  tp[N}[N}TdA)[u]  +  {Su}T{F}         (2.15) 

J  a  J  A 

By  invoking  the  principle  of  virtual  displacements  and  noting  that  {Su}  is 
arbitrary,  we  get  the  equilibrium  equations  in  the  following  form; 


[K}{u}  +  [M]{u}  =  [F]  +  [F}d  (2.16) 


where. 


[M\  =  I  tp[N)[N}TdA  (2.17) 

J  A 

ft   1 


[F\  = 


Fn 


(2.18) 


[F]d  =  J[N]Tds  (2.19) 

It  may  be  noted  that  the  matrix  [M]  is  the  mass  matrix,  [F]  is  the  vector  of 
point  loads  and  [F]    is  the  vector  of  consistent  loads. 


C.     LAMINATE  THEORY 

1.  Introduction 

In  the  expression  for  [K]  matrix,  there  are  three  unknown  matrices.  In 
this  section,  we  discuss  about  calculation  of  [Q]  matrix.  The  matrix  [Q]  ,  relat- 
ing stresses  and  strains,  consists  of  material  stiffness  coefficients.  It  reflects  the 
properties  of  both  fibers  and  matrix. 

The  behavior  of  laminated  plates  is  characterized  by  possible  coupling 
between  membrane  action  and  bending  action. 

Following  discussions  review  some  aspects  of  the  laminate  theory. 

2.  Strain-displacement  Relations 

As  shown  in  Figure  2.1,  a  general  strain  field  converts  configuration  012 


to  0T2 


Linearized  normal  strains  may  be  written  as 

du 


dx  <2-20> 


and 


*  =  |  (2.2!) 

The  shear  strain  is  defined  as  the  amount  of  change  in  a  right  angle.  For 
small  angles, 


7*y    =    /?i+/?2  (2.22) 

du      dv 
=     Ty+T*  (2-23) 

Similar  expressions  may  be  derived  for  other  strain  components. 


The  foregoing  strain-displacement  relations  can  be  summarized  in  a  matrix- 
operator  form. 
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3.     Stress-Strain  Relations 

The  generalized  Hooke's  Law  takes  the  form, 


M  =  IQM 


where,  the  stress  vector  is  given  by 


(2.24) 


(2.25) 


{a}  =  {ax  oy  az  rxy  ryz  rzx\ 
and  the  strain  vector  is  given  by 


(2.26) 


{e}  =  {ex  ey  ez  ^xy  -yyz  ^zx]  (2.27) 

The  material  stiffness  matrix  [Q]  ,  contains  nine  independent  coefficients. 
An  orthotropic  material  has  only  five  independent  elastic  coefficients  —Ex,E2,  V\2i  U2i,and  G\2- 
More  explicitly,  Equation  (2.25)  may  be  written  as 


Txy 

Tyz 
TZx 


>    = 


Gii  Q12  Qis  o 

Qn  Q22  Q23  0 

Q31  Q32  Q33  0 

0  0  0  Q44 

0  0  0  0 

0  0  0  0 


0 
0 
0 
0 

£55 

0 


/                    * 

0 

£x 

0 

ev 

0 

< 

tz 

0 

Ixy 

0 

lyz 

Qee 

k  Izx    t 

(2.28) 


*    u+  f d^ 


»3 


*} 


* L_ 


0 


e-     14      — * 


<- 


dx 


2' 


9l    w-Sdat  k 


V  +  fdx 


Figure  2.1:   Displacement  and  Distortion  of  Differential  Lengths  dx  and 
dy. 
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For  planar  orthrotropic  material,  the  stress-strain  relation  reduces  to: 


where, 


°y 

TXy 

Ty* 


>  = 


Qii  Qu  0  0  0 

Q21  Q22  0  0  0 

0  0  Q44  0  0 

0  0  0  Q55  0 

0  0  0  0  Qee 


tx 

ty 

Ixy 
7yz 


{    lz 


(2.29) 


Qu 

— 

1  -  ^12  ^21 

Q22 

= 

E2 

1  -  UU  U2i 

Q44 

= 

G\2 

Qss 

= 

G23 

Qee 

z=. 

G13 

(2.30) 

(2.31) 
(2.32) 
(2.33) 
(2.34) 

It  may  be  noted  that  shears  ryz  and  rzx  are  obtained  to  account  for 
transverse  shears. 

4.     Lamina  of  Composite  Materials 

Composite  structures  are  built  of  individual  lamina,  which  are  stacked 
into  several  number  of  layers  to  form  a  laminate.  Each  lamina  consists  of,  typically, 
uniaxial  fibers  embedded  in  a  matrix,  such  as  a  resin.  In  Figure  2.2,  the  principal 
material  axes  are  labelled  1  and  2,  that  is,  1-direction  is  parallel  to  the  fibers  direction 
and  the  2-direction  is  normal  to  them. 

It  may  be  noted  that  in  each  lamina  there  exists  a  state  of  plane  stress. 

The  state  of  stress  is  also  shown  in  the  Figure  2.2,  in  both  1-2  and  x-y 
coordinate  system.  The  computation  stresses  in  different  coordinate  system  follows 
the  usual  transformation  rules  [Ref.    21].    We  follow  the  notation  that  01,03  are 
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6+ 

X 


<5+ 

X 


*y+ 


Figure  2.2:  Lamina  Coordinate  System  (2-Dimension) 
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normal  stresses,  03,04  and  05  are  the  shear  stress  in  1-2  system.  ei,e2,  63,64  and  e5 
are  the  corresponding  strains  in  1-2  system. 


=  Plm 

where  the  transformation  matrix  is  given  by 


02 
<r4 


'xy   . 


\T\cl  = 


m  n        2mn 

™2     -2mn 

,2 


n  m 

[  -mn    mn     (m2  -  n2)  J 

The  direction  cosines  of  the  unit  normal  are  determined  from 


(2.35) 


(2.36) 


m 


=  COS0 


(2.37) 


n  =  sin# 


(2.38) 


The  subscript  CL  refers  to  two-dimensional  case,  that  is  the  x-y  plane 
only.  Similarly,  the  strains  are  related  in  the  two  coordinate  systems  by 


(2.39) 


By  including  transverse  shears,  we  modify  these  transformations  as  fol- 
lows: 


€l 

tx 

^2 

=  \T\cl 

ey 

£12   . 

.  7xi/  . 

<J\ 

Gx 

°2 

°y 

<*3 

=  [T) 

TXy 

04 

TVZ 

.  a*  . 

'zx 

(2.40) 
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t\ 

€x 

e2 

% 

e3 

=  [T] 

Ixy 

e4 

lyz 

.   65   . 

"Yzx 

(2.41) 


where  the  transformation  is  given  by, 


[71  = 


m2 

n2 

2mn 

0 

0 

n2 

m2 

—Imn 

0 

0 

—mn 

mn 

2           2 

m   —  n 

0 

0 

0 

0 

0 

m 

— n 

0 

0 

0 

n 

m 

(2.42) 


The  stresses  and  strains  in  x-y  system  may  simply  be  obtained  from  1-2 
system  by  inverting  the  matrix  [T]. 


ax 

(             \ 

°y 

Txy 

=  IT]'1  < 

a2 

&3 

Tyz 

a4 

>zx 

>  °s  t 

(2.43) 


and, 


tx 

Cl 

£y 

Ixy 

=  PI'1 

C2 
C3 

lyz 

e4 

i ' zx 

.   65   . 

(2.44) 


with 


pr  = 


m2 

n2 

— 2ran 

0 

0 

n2 

m2 

2mn 

0 

0 

mn 

—mn 

m2  —  n2 

0 

0 

0 

0 

0 

m 

n 

0 

0 

0 

— n 

rr 

(2.45) 


The  stress-strain  relations,  then,  in  x-y  system  assumes  the  following 


form, 
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<jx 

Qn 

Ql2     Ql3     0 

0 

tx 

°y 

Q21 

Q22  O23  0 

0 

ev 

TXy 

►   = 

Q31 

Q32  Q33  0 

0 

< 

Ixy 

Tyz 

0 

0        0        Q44 

0 

lyz 

Tzx     . 

0 

0        0        0 

Q55 

„    Izx 

(2.46) 


where, 


1-1 


=  PT»]P1 


(2.47) 


5.      Laminate  Analysis 


As  discussed  earlier,  the  laminae  are  stacked  to  obtain  a  laminate.  In  this 
section,  the  theory  associated  with  the  mechanics  of  laminates  is  described. 

Consider  a  laminate  composed  of  TV  lamina.  For  the  kth  lamina  of  the 
laminate,  Equation  (2.46)  can  be  written  as  follows: 


<7„ 


'  xy 


'yz 


-Ml. 


K 


Ixy 
lyz 
L    Izx 


(2.48) 


A' 


where  all  matrices  must  have  the  subscript  K  due  to  orientation  of  the  particular 
lamina  with  respect  to  the  plate  x-y  coordinates. 

The  functional  form  of  the  displacement  for  a  plate  are  given  by: 


u(x,y,z)    =    u0(x,y)  +  za(x,y) 
v(x,y,z)    =    v0(x,y)  +  z(3(x,y) 
w(x,y)    =    w0{x,y) 


(2-49) 
(2.50) 
(2.51) 


where  u0,  v0  and  w0  are  middle  surface  in  plane  displacements,  a  and  /?  are  related 
to  the  rotations. 
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Substituting  Equation  (2.49),  (2.50)  and  (2.51)  into  the  Equation  (2.24) 
results  in: 


dv0        dp 

«•  =  i£  +  %  <2-53> 

ez    =    0  (2.54) 

1  (du0      dv0\ 

€-    =    2^  +  ^J  (2-55) 

*    m    2{0+Ty)  (256) 

The  mid  surface  strains  are  given  by  the  following  relations: 

V.    =    £  (2.58) 

V.    =    ^  (2-59) 

<~  =  2^+a7J  (2'60) 

The  curvature  terms,  associated  with  transverse  bending  are  written  in 
the  form, 


K,    =    |  (2.6!) 

«y    =    7T  (2-62) 


dy 

1  /da      ojT 

2  \%  +  fo 


(2.63) 


On  substituting  the  strain-displacement  relations  into  the  stress-strain 
relations,  we  obtain  stresses  in  terms  of  displacement  components, 
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<yx 

txo 

+ 

ZKX 

°y 

£Vo 

+ 

ZKy 

TXy 

>     = 

[Q\ 

k\ 

lxy0 

+ 

ZKXy 

Tyx 

lyz 

'zx    J 

K 

„    Izx 

(2.64) 


A' 


For  plate/shell  type  structures,  we  define  stress  resultants  Nx,  Ny,  Nxy, 
Qx,  Qy,  Mx,  My,  and  Mxy,  as  shown  in  Figure  2.3.  For  the  kth  layer,  we  may  write, 


f  N*  } 

°x 

Ny 

■*■  *  xy 

Qx 

ft/2 
>   =    /           < 
J-t/2 

°y 

<TXy 

°yz 

I  Qv    \ 

b    &zx    t 

dz 


(2.65) 


For  a  laminate  composed  of  N-layers,  the  normal  stress  resultants  are 


given  by 


Nx 
N„ 


N. 


xy 


N      ttk 


xy 


dz 


or,  in  terms  of  strains,  we  have 


Nx 
N. 


xy 


N  .tk  x°  rtk 


Kx 

i 

Ky 

zdz  y 

.   Kxy   _ 

, 

The  equation  (2.67)  may  be  written  in  the  form 


(2.66) 


(2.67) 


[N]  =  [A]  [e0]  +  [B]  M 
with  the  elements  of  the  matrices  [A]  and  [B]  given  by 


(2.68) 


N 
k=l 

zfc=l 


(2.69) 
(2.70) 
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The  moment  resultants  are  obtained  for  a  K     layer  as  follows 


'  Mx 


M. 


xy 


-r. 


t/2 
t/2 


°x 

&xy 


►  zdz 


(2.71) 


or, 


[M]  =  [B]  [e0]  +  [D]  [k] 

1  N 
where  DXJ     =  -  £  (<50)^  (^  -  ^.^ 


(2.72) 
(2.73) 


k=i 


D.     [JB],  THE  STRAIN-DISPLACEMENT  MATRIX 

1.  Introduction 

In  this  section,  we  describe  how  the  matrix  [B]  may  be  calculated. 
The  matrix  [B],  which  relates  the  strains  and  displacements  at  nodal 
points  of  an  element  may  be  obtained  in  the  form, 

{e}  =  [B]{u}  (2.74) 

It  may  be  observed  that  this  matrix  depends  on  the  choice  of  the  nodal 
degrees  of  freedom,  the  shape  function  and  the  form  of  strain-displacement  relations. 

In  the  discussion  that  follows,  we  describe  the  nine  noded  bi-quadratic 
Lagrangian  isoparametric  element.  There  are  five  degrees  of  freedom  associated  with 
each  node. 

2.  Shape  Functions  and  Their  Derivatives 

A  finite  element  is  called  an  isoparametric  element  if  the  same  interpola- 
tion functions  define  both  the  geometry  and  displacements. 

Figure  2.4  shows  the  element  in  the  mapped  space.  For  the  in-plane 
deformation,  the  geometry  may  be  interpolated 
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Figure  2.3:  Positive  Directions  for  Stress  Result 
tor  a  Lamina 


ants  and  Stress  Couples 
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y 


=  [N]{x1y1x2y2  ••■x9y9} 


(2.75) 


and  the  displacements  given  by 


=  [TV]  {Ux  VXU2V2  •  •  •  U9  **>} 


(2.76) 


where  the  matrix  of  shape  functions  is  given  by 


[N)  = 


JVi    0      N2    0       ...    N9    0 
0       Nx    0      N2    •••    0       7V9 


(2.77) 


We  may  also  write  the  displacements  in  a  summation  notation  as 


u 

i=l 

V 

=  Eiv,^, 

(2.78) 
(2.79) 


t=i 


In  the  case  of  plate  bending,  three  more  degrees  of  freedom  are  added  to 
the  planar  displacements. 

The  transverse  displacement  w  ,  9X  and  0y  are  the  rotations  of  the  normal 
to  the  undeformed  middle  surface  in  the  x  —  z  and  y  —  z  plane,  respectively. 


(2.80) 
(2.81) 
(2.82) 


w 

0X 

Oy 

=  E% 

«=1 


The  shape  functions  and  their  derivatives,  which  are  needed  for  the  com- 
putations of  strains,  are  given  in  Table  2.1. 
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Figure  2.4:  Lagrangian  Plate  Element 
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Figure  2.5:  Degrees  of  Freedom  at  a  Typical  Node  i. 

It  may  be  noted  that  the  rotational  degrees  of  freedom  are  treated  as 
independent  quantities,  following  Mindlin  theory  [Refs.  37,  38,  39]. 
3.      Jacobian  Transformation  Matrix 

The  chain  rules  for  differentiation  of  A'(£,  rj)  with  respect  to  £  and  77  gives, 


dNj      ON,  dx      dNjdy 


(2.83) 


or,  in  matrix  form, 


dNj      dNjdx      ON,  dy 

dr]         dx  di]       dy  dr) 


(2.84) 


dx     dy 

dx     dy 
di)     dij 
The  Cartesian  derivatives  are  given  by 


"  dNi  ' 

dNi 

dr] 

= 

r  dNi 

L  dy 


(2.85) 
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Figure  2.6:  Displacements  in  the  x-z  Plane 
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TABLE  2.1:  Shape  Function  and  Their  Derivatives 
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dNj 
fa 


dx     dy  1  -1  f  dNj 

MM  di 

dx     dy  dNi 


(2.86) 


dy    J       [  di)     drj  ]       [    dr\ 
The  Jacobian  matrix  [J],  which  contains  the  derivates  of  the  global  coor- 
dinates with  respect  to  the  local  coordinates  is  given  by, 


in- 


dx     dy 
dx     dy 


(2.87) 


drj     drj 

By  using  isoparametric  element  concepts  [Ref.   1],  we  may  write  the  Ja- 
cobian elements  as, 


7  f  8N> 


(2.88) 


T  r9Ni 


(2.89) 


j  -r™ix 
21  ~fe&?  ' 


(2.90) 


J22  =  2^  -Q-y»- 

The  inverse  of  the  Jacobian  matrix  can  be  expressed  as 


(2.91) 


[J]-1  = 


1 


J 


J22 

— J21 


-J12 

J11 


with  the  determinant  of  the  Jacobian  given  by 


(2.92) 


J  \—  J11J22  ~  J21J12 


(2.93) 
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The  Cartesian  derivatives  are,  then,  given  by 


dNj 

dx 


I  j  I  v    di 


7   dNi       _   dNt 


dr\ 


(2.94) 


dNi 


1 


4. 


Formation  of  [B]  Matrix 


—  ^21-^T"  +  «/ll' 


ac 


$77 


The  equation  (2.74)  may  be  written  as 


'  ex 

r  — 

3x 

0    1 

«» 

= 

0 

d 

* 

.    *«v  . 

L  av 

3x    - 

(2.95) 


(2.96) 


Replacing  the  displacements  in  terms  of  nodal  displacements,  we  get 


or,  in  matrix  form 


cy 

Ixy 


•  Tfi    0 

0  ^ 

dy  dx 


Vi 


(2.97) 


{<}  -  [5,1  {"} 


(2.98) 


where  [B{]  is  given  by 


[Bi]  = 


For  plate  bending,  we  have, 


3N, 
9x 

0 

dN, 
dy 


0 

3AT. 

37^, 
dx 


(2.99) 


ey 

7xy 

or,  in  terms  of  nodal  displacements, 


0    fx  o 

0    °  h 

o  i  i 

3y  3x 


If) 


(2.100) 
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The  transverse  shears  are  expressed  by 


tx 

0 

dx 

€v 

= 

0 

0 

Ixy   . 

0 

dy(N, 

&xi 


(2.101) 


or, 


lyz 
~1zx 


dy 

f    1    o 

ox 


(2.102) 


~izx 


3x 


0       -EiV, 
-EAi       0 

We  may  write  [B{]  for  plate  elements  as 


(2.103) 


[Bi]  = 


0 
0 
0 

dNj 
dy 

dNi 

dx 


^Vl       o 
dx 


0 


dNj 

dy 


0 

-Ni 


dNj 

dftj 
dx 

-Ni 


(2.104) 


E.     GAUSSIAN  QUADRATURE 

1.  Introduction 

In  this  section,  we  discuss  aspects  of  numerical  integration.  Here  we 
describe  the  Gaussian  method  to  calculate  the  integrals,  which  is  widely  used  in 
finite  element  work  [Refs.  38,  39]. 

2.  Summary  of  Gauss  Quadrature 
To  approximate  an  integral  J, 

1=1    <f>di  (2.105) 

where  <f>(()  is  a  function  defined  in  that  region.    As  shown  in  the  Figure  2.7,  we 
sample  <f>  at  the  midpoint  of  the  interval  and  multiply  by  the  length  of  the  interval. 


27 


Generalization  of  Equation  (2.105)  leads  to  the  formula 


/  =  E^(W  (2-106) 


/  =  Wi^fi)  +  Witfjh)  +  •  •  •  +  WmtfjU)  (2.107) 

where  £,  is  the  location  of  the  integration  point  i  with  respect  to  the  origin,  W;  is 
the  weighting  factor  for  point  i  and  n  is  the  number  of  points  at  which  <^>(£)  is  to  be 
evaluated. 

The  sampling  points  and  weights  for  Gaussian  Quadrature,  which  is 
adapted  in  the  present  study,  are  listed  in  Table  2.2. 

In  two  dimensions,  we  may  approximate  /  by 

/  =  EEW(t.1j)  (2-108) 

«'       3 

In  Equation  (2.6)  and  (2.7)  each  coefficient  in  the  integrand  matrix  must 
be  integrated  as  described  above. 
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Figure  2.7:    Gaussian  Quadrature  (£  Coordinate)  using  Three  Sampling 
Points 
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n 

±f. 

R, 

I 

00 

20 

1 

0  5773502692 

1  0 

3 

0  7745*^-692 
0  0 

0  5555555556 
0  88888XS8X9 

4 

0  K61I '63116 
0  33<W81»)436 

0  3478548451 
0  6521451549 

5 

0  9itf.  1798459 
0.53M693I0I 
00 

0  236V26X85I 
0  47x62X6705 
0.56XXXX8889 

6 

0  9324MJSI42 
0  66!2(>9'n65 
0.2386191861 

0.1713244924 
0  ^607615730 
0  4679139346 

7 

0  949KP9I23 
0  74 1531 1856 
0  405845I5M 

0  (» 

0  I2'U849662 
0  2797053915 
0  38 18 '00505 
0  4 17959 1837 

8 

0  9602X98565 
0  "Nt>6664774 
0  525V24099 
0  |X34U(>425 

0.1012285363 
0  2223810345 
0  3M706M59 
0  3626837834 

TABLE  2.2:  Coefficients  for  Gaussian  Quadratu 


re 
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III.  PROGRAM  IMPLEMENTATION 

A.     INTRODUCTION 

After  having  formed  the  element  stiffness  matrices,  the  global  matrices  are 
formed  in  a  standard  manner  [Refs.  1,  38],  resulting  in,  for  static  analysis, 


where     {F} 

W 
{u} 


{F}  =  [K]  {u} 

External  loads  vector 

Stiffness  matrix 

Nodal  displacements  vector 


If  we  know  {F}  and  [K],  {u}  may  be  computed  for  a  static  problem.  Typical 
steps  for  static  analysis  may  be  described  as  follows: 

Step  1:  Input  the  material  properties,  plate  coordinates,  loads,  the  number  of 
integration  points,  boundary  conditions. 

Step  2:  For  the  composite,  input  number  of  layers,  elements  and  fiber  orien- 
tations. 

Step  3:  Determine  element  matrices  in  global  coordinate  system. 

Step  4:  Assemble  the  element  matrices. 

Step  5:  Solve  for  displacements  and  stresses. 

The  Figure  3.1  shows  the  flow  chart  of  the  program. 

B.     SOLUTION  PROCEDURE 

In  order  to  obtain  numerical  solution,  the  element  matrices  were  programmed 
and  incorporated  into  an  existing  finite  element  program.  The  implementation  in 
double  precision  was  done  on  a  32-bit  Apollo  D-3000  series  computer. 

Following  steps  describe  the  element  subroutine. 
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BEGIN 


INPUT/ 
N=  1, NO. OF  LAYERS 


CALCULATE       Q 


/ 


INTEGRATION  (2*2)  OR  (3*3) 


SHAPE  FUNCTION 
JACOBIAN  MATRIX 


B    MATRIX 
FORMULATION  OF       K 


K 
T 


END 


Figure  3.1:  Flow  Chart  (Main  Program) 
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Step  1:  Material  properties  from  main  program  are  read. 

Step  2:  Calculate  [Q]  according  to  number  of  layers. 

Step  3:  Select  integration  point  2x2  or  3x2  or  3x3. 

Step  4:  Establish  shape  function. 

Step  5:  Determine  the  Jacobian  matrix  and  [B]. 

Step  6:  Formulation  of  [K]. 

Step  7:  For  each  integration  point,  do  steps  4,  5,  and  6,  and  accumulate  [K], 

Step  8:  Calculate  [K]  for  each  element. 

Step  9:  Return  to  main  program. 

The  flow  chart  for  element  level  computation  is  given  in  Figure  3.2. 
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INPUT. 
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CALCULATE       Q 
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INTEGRATION  (2*2)  OR  (3*3) 


SHAPE  FUNCTION 
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B    MATRIX 
FORMULATION  OF       K 


JL 


K 


END 


Figure  3.2:  Flow  Chart  (KQML9) 
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IV.  NUMERICAL  EXAMPLES 

A.  INTRODUCTION 

This  section  discusses  the  validation  of  the  element  formulation  and  imple- 
mentation by  means  of  selected  numerical  examples.  The  isotropic  elastic  plates  are 
solved  under  different  boundary  conditions  and  loadings. 

This  is  followed  by  application  to  selected  laminated  composite  plates  to  check 
the  formulation  for  such  applications.  The  effects  of  thickness  and  integration 
schemes  also  are  investigated. 

B.  SIMPLY  SUPPORTED  SQUARE  ISOTROPIC  PLATE 

1.      Simply  supported  plate  under  concentrated  load. 

Consider  a  rectangular  isotropic  plate  under  a  central  point  load.  The 
geometry  and  the  boundary  conditions  are  shown  in  the  Figure  4.1.  The  structure 
is  modeled  using  the  double  symmetry. 

The  properties  of  material  /  are  given  by 

E     =     10.92  xlO6  (psi) 

v     =     0.3 

L     =     10  in 

P     =     400  (lbs) 
The  geometric  boundary  conditions  of  the  simply  supported  plate  are 

w  =  0  on  all  edges. 

As  shown  in  the  Figure  4.2,  results  depicting  maximum  displacement 

obtained  using  2x2  integration  points,  3x3  integration  points  and  'heterosis'  [Ref. 

4]  elements  are  compared.  As  the  thickness  decreases,  the  element  is  more  effective 

for  integration  schemes  and  for  thick  plates,  the  error  is  approximately  20%.   It  is 

possible  that  by  increasing  the  number  of  elements,  it  may  improve  the  solution. 
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Figure  4.1:    A   Rectangular  Isotropic  Plate  under  Concentrated   Point 
Load 
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WlWanai  vs.  No.  of  elements  are  shown  in  Figure  4.3.  As  the  number  of 
elements  are  increased,  the  solution  converges  to  the  analytically  predicted  values. 
It  may  be  noted  that  reduced  order  of  integration  is  more  flexible,  and  approaches 
the  analytical  solution  as  an  upper  bound. 

2.     Simply  Supported  Plate  under  Distributed  Load 

Next  we  consider  a  rectangular  plate  under  a  distributed  load.  The  ge- 
ometry and  the  boundary  conditions  are  the  same  as  in  the  previous  example.  The 
model  of  the  structure,  using  the  symmetry,  is  shown  in  Figure  4.4. 

The  material  properties  for  this  example  are  given  earlier  as  material  I. 

The  normalized  maximum  deflection  is  shown  plotted  against  the  number 
of  elements  in  Figure  4.5.  It  may  be  observed  that  different  integration  schemes 
converge  towards  a  lower  bound.  Figure  4.6  depicts  the  effectiveness  of  this  element 
in  predicting  behavior  of  thick  and  thin  plates. 

This  element  gives  better  predictions  for  thick  plates  than  heterosis  ele- 
ment, however,  for  thin  plates,  heterosis  element  appears  to  be  better. 

C.     CLAMPED-CLAMPED  SQUARE  ISOTROPIC  PLATE 

Consider  a  rectangular  plate  clamped  on  all  sides,  under  a  distributed  load. 
The  geometry  and  boundary  conditions  are  shown  in  Figure  4.7.  The  boundary 
conditions  for  this  problem  are  implemented  by  prescribing  all  degrees  of  freedom 
to  zero,  on  all  four  edges.  The  structure,  again,  is  modeled  using  the  symmetry. 

The  material  properties  for  this  plate  are  same  as  for  material  I. 

The  results  showing  the  normalized  maximum  displacements  are  shown  in 
Figure  4.8.  The  heterosis  element  results  in  about  ten  percent  error  while  the 
Lagrangian  element  gives  about  twenty  to  twenty-three  percent  for  thick  plates. 
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Figure  4.2:    \V/Wanal  vs.    Log  (L/T)  on  Simply  Supported  Square  Plate 
under  Central  Point  Load. 
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Figure  4.3:    H7Wana/  vs.    No.    of  Elements  on  Simply  SupporteoTSquare 
Plate  under  Central  Point  Load. 


39 


ss 


ss 


Figure  4.4:  A  Rectangular  Isotropic  Plate  under  Distributed  Load 
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Figure  4.5:    W/W(anai)  vs-    No.    of  Elements  on  Simply  Supported  Plate 
Under  Distributed  Load 
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Figure  4.6:    W/W(anal)  vs.     Log{L/T)  on  Simply  Supported  Plate  Under 
Distributed  Load 
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Figure  4.7:    A  Clamped-Clamped  Rectangular  Plate  Under  Distributed 
Load 
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However,  for  thin  plates,  the  predictions  are  about  the  same  using  either  elements. 
Figure  4.9  shows  the  convergence  characteristics  of  the  element. 

D.  CANTILEVERED  ISOTROPIC  PLATE 

Next  example  considered  is  a  cantilevered  plate.  The  geometry  and  boundary 
conditions  are  shown  in  the  Figure  4.10.  The  material  properties,  (material  II),  are 
given  by 

E  =  1  x  I04(k/in2)  (4.1) 

v  =  0.3  (4.2) 

L  =  lO(in)  (4.3) 

t  =  0.2(m)  (4.4) 

The  result  of  this  example  is  shown  in  the  Figure  4.11.  As  the  number  of 
elements  is  increased,  the  results  converge,  for  all  integration  schemes. 

E.  SIMPLY  SUPPORTED  LAMINATED  PLATE 

1.      Graphite-epoxy 

Consider  a  rectangular  composite  plate  under  a  sinusoidal  load.    Two 

types  of  construction  are  used. 

Case  I      :     Graphite- Epoxy     0o/90°/90o/0° 
Case  II     :     Graphite- Epoxy     0°/  -  60°/60°/0 
Material  properties  are  given  by 


Gu     =    0.6  (4.5) 


o 


E2 

=    0.5  (4.6) 


G22 


E2 

vu    =    0.25  (4.7) 
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Figure  4.8:    W/W{ana[)  vs.     Log(L/T)  for  Clamped-Clamped  Plate  under 
Distributed  Load 
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Figure    4.9:       Nondimensional    Deflection    vs.       No.       of   Element     for 
Clamped-Clamped  Plate  under  Distributed  Load 
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Figure  4.10:  Cantilevered  Plate 
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Figure  4.11:  VK/VV(ana0  vs.  No.  of  Element  for  Cantilevered  Plate 
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f-    =    40  (4.8) 

The  plate  is  subjected  to  a  transverse  sinusoidal  load  of  intensity 


q  =  q0  sin  —  sin  —  (4.9) 

This  problem  is  studied  to  compare  finite  element  results  with  those  of 
shear  deformation  theory  [Ref.  24]  and  elasticity  solution  [Refs.  16,  17,  18].  The 
results  are  in  Figure  4.12  for  Case  I.  The  deflection  obtained  from  the  finite  element 
solutions  agree  very  well  with  the  exact  elasticity  solutions.  The  solutions  agree 
well  for  thin  plates,  while  a  large  discrepancy  between  the  present  predictions  and 
Reference  24  is  observed. 

The  maximum  deflection  is  plotted  as  the  number  of  elements  are  in- 
creased in  Figure  4.13.    The  convergence  trend  may  be  noted  as  the  number  of 
elements  are  increased.  In  Figure  4.14,  the  results  for  the  second  type  of  composite 
(0°/  -  6076070°)  is  depicted. 
2.      Glass-Epoxy 

This  section  considers  rectangular  composite  plate  under  sinusoidal  load. 
The  two  types  of  construction  are  Case  III,  0°/90o/90°/0  and  Case  IV,  0°  /  - 
60°/60°/0°.  The  plates  are  square  and  simply  supported. 

The  material  properties  are  given  by 


if  =  °-6  <4-10) 

IT  =  0.5  (4.11) 

1/12  =  0.25  (4.12) 

#-  =  3  (4.13) 
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Figure  4.12:  Central  Deflection  Wmax  of  4- Layered  Square  Plate 

(Case  I:  Graphite-Epoxy) 
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Figure  4.13:  Wmax  vs.  No.  of  Elements  for  4-Layered  Square  Plate 

(Case  I:  Graphite-Epoxy) 


51 


• 

If 

o 

t 

j 

i 

001 


OS  0  9 

„  0I«(XVW)ii 


N 


z 

w 

w 

o 
o  6 


OM 


Figure  4.14:  lVmox  vs.  No.  of  Elements  for  4-Layered  Square  Plate 

(Case  II:  Graphite- Epoxy) 
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The  plate  is  subjected  to  a  transverse  sinusoidal  load  of  intensity. 

.    ttx    .    Try 
q  =  q0sm—sm—  (4.14) 

Figure  4.15  shows  the  maximum  deflection  of  this  element  and  of  shear 
deformation  theory  [Ref.  24]  and  elasticity  solution  [Refs.  16,  17,  18]  for  analytical 
solution.  Good  agreement  is  obtained  for  thin  plates  while  discrepancies  exist  for 
thick  plates.  Convergence  characteristics  are  shown  in  Figures  4.16  and  4.17. 
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Figure  4.15:  Central  Deflection  Wmax  for  4-Layered  Square  Plate 

(Case  III:  Glass-Epoxy) 
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Figure  4.16:  Wmax  vs.  No.  of  Elements  for  4-Layered  Square  Plate 

(Case  III:  Glass-Epoxy) 
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Figure  4.17:  VVmax  vs.  No.  of  Elements  for  4-Layered  Square  Plate 

(Case  IV:  Glass-Epoxy) 
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V.  CONCLUSIONS 

This  study  is  directed  towards  understanding  the  linear  static  analysis  of  plates 
composed  of  both  isotropic  and  laminated  composites. 

The  formulation  is  based  on  the  principle  of  virtual  work.  A  finite  element 
formulation  is  presented  as  a  model  for  the  analysis  of  laminated  anisotropic  plate 
structures.  Several  numerical  examples  for  the  isotropic  elastic  plates  are  solved  for 
different  boundary  conditions  and  loadings. 

The  results  show  that  bi-quadratic  isoparametric  lagrange  plate  bending  ele- 
ment is  effective  for  analysis  of  laminated  anisotropic  plates.  Numerical  solutions 
agree  well  with  analytical  solution.  Further,  it  is  observed  that  as  the  number  of 
elements  are  increased,  the  convergence  to  analytical  solution  is  assured. 

The  element  predicts  good  results  for  thin  plates  while  large  discrepancies  exist 
for  thick  plates  and  calls  for  further  investigation.  In  general,  a  3  x  3  integration 
seems  to  predict  the  solutions  well. 

More  numerical  experiments  need  to  be  done  regarding  selective  integrations, 
and  apply  to  a  variety  of  layer  configurations  for  composite  plates.  Another  area  is 
in  including  the  nonlinear  terms  for  buckling  problems.  The  free-vibration  analysis 
to  predict  the  natural  frequencies  and  mode  shapes  is  an  obvious  extension. 
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